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CONTRACTOR REPORT 


A FINITE ELEMENT COMPUTATIONAL METHOD FOR 
HIGH REYNOLDS NUMBER LAMINAR FLOWS 

I. INTRODUCTION 


The finite element computational methods for the Navier- Stokes equations using 
the primative variables of velocity and pressure can, in general, be categorized into 
three groups. These are the velocity-pressure integrated mixed interpolation methods, 
the penalty methods, and the segregated velocity-pressure solution methods. 

In the velocity-pressure integrated mixed interpolation methods, the order of 
interpolating polynomial for velocity is chosen to be one order higher than that of 
pressure. For example, if constant elements are used for pressure, then linear 
elements are used for velocity; and if linear elements are used for pressure, then 
quadratic elements are used for velocity. Unfortunately, most of the existing velocity- 
pressure integrated, mixed interpolation methods exhibit spurious pressure modes 
which become more severe as the Reynolds number is increased, and convergent 
solution can not be obtained for high Reynolds number flows. Details of the method 
can be found in Taylor and Hughes [1] and the references therein. 

In the penalty method, the pressure variable is pre- eliminated from the Navier- 
Stokes equations by penalizing the conservation of mass equation, and hence the con- 
tinuity condition is satisfied approximately. If necessary, pressure can be recovered 
using the penalized conservation of mass equation in the post process. Unfortunately, 
most of the computational results obtained using the penalty methods revealed that 
the conservation of mass equation is not satisfied rigorously so that the computational 
results deteriorate as the Reynolds number is increased. Details of various penalty 
methods and the computational results can be found in Zienkiewicz et al. [2], Engelman 
et al. [3], Kikuchi et al. [4], and Heinrich and Marshall [5], among many others. 

Due to the shortcomings of the previous two classes of methods, and partly 
influenced by the success of the finite difference computational methods of segregated 
formulation of the N avier- Stokes equations, such as the SIMPLE (Semi- Implicit Method 
for Pressure-Linked Equations) algorithm [6], the finite element versions of segregated 
formulation of the Navier- Stokes equations began to appear in recent years [7,8,9]. 

But the computational results thus obtained did not show any significant improvement 
over the previous two classes of methods. 

The finite element method has certain advantages over the finite difference 
methods. These advantages are: capability to model complicated domain more pre- 

cisely, capability to include various boundary conditions more naturally, and capability 
to show the convergence nature of the method mathematically. For low Reynolds 
number flows, a number of example cases can be found for which the finite element 
methods yielded more accurate results than the finite difference methods [ 10] . But 
for high Reynolds number flows, the finite difference methods yielded more accurate 
results, i.e., the finite difference computational results compared more favorably with 
experimental data than those of the finite element methods, especially when pressure 
driven recirculation zones exist in the flow field. Thus, the advantages of the finite 
element method have not been well demonstrated for high Reynolds number viscous 
flows as yet. 



In this report, a new velocity- pressure integrated, mixed interpolation Galerkin 
finite element method for the Navier- Stokes equations is presented. The finite element 
system of equations is obtained using the Galerkin finite element method, and the 
system of equations was solved using a frontal solver [1,11]. The present method 
yielded accurate computational results for low Reynolds number flows as well as high 
Reynolds number flows. It was also found that the method can capture subtle pres- 
sure driven recirculation zones, and that the method did not yield spurious pressure 
modes . 


II. FINITE ELEMENT EQUATIONS 


A finite element computational method for two- and three-dimensional laminar, 
steady , incompressible flows is described below . The method is based on the inte- 
grated velocity-pressure formulation of the Navier- Stokes equations using a new mixed 
interpolation of velocity and pressure. 

In the following discussions, consistent notations are used throughout, and 
repeated indices imply summation over the indices, unless otherwise specified. 


2-1. The Navier- Stokes Equations 

The form of the Navier- Stokes equations used herein are given as: 
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ft is the open bounded domain of the problem, the subscripts i and j denote coordinate 
directions, p is the density of the fluid, u. is the velocity component in the i-th 

coordinate direction, p is the pressure, y is the molecular viscosity of the fluid, b. 

is the body force in the i-th coordinate direction, t.. is the stress tensor, e.j is the 
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strain rate tensor, and Sy is the Kronecker delta such that 6y. = 1 for i = j and 
6y = 0 for i ^ j. The boundary conditions used are given as: 

u = u (x) for x e 3 ft, (5) 

T. = t-.. n. for x e 3 fi, (6) 

1 1J J L 

3^2 = 9 

(7) 

3 ^ n 9 £7 2 = $ 

where x = (x,y) for 2-D case and x = (x,y,z) for 3-D case, 9 ft is the boundary of 
the domain, is part of the boundary on which Dirichlet boundary condition is 

specified, 9^ 2 is another part of the boundary on which Neumann boundary condition 

is specified, is the surface traction, and <|> is the null space. A Dirichlet pressure 

boundary condition is specified at an arbitrary pressure node in the flow domain. 

2-2. Method of Weighted Residuals and the Galerkin Finite Element Method 

In the context of the method of weighted residuals, the test function for the 
momentum equation and the conservation of mass equation are denoted as W u (x) and 

W (x), respectively. The weak form of the Navier-Stokes equations are obtained to be 

p ~ 



Integrating by parts of the term containing the stress tensor in equation (8) 
yields : 
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where l W T..n.ds term has been dropped out since the test function W (x) vanishes 
J u 1 ] j u ~ 

9 

on the boundary 3 

In order to apply the finite element method to the weak form of the Navier- 
Stokes equations, equations (9) and (10), the global domain is discretized into a 
number of closed bounded subregions, called the finite elements (ft g ), such that 


l = !1 J 32, 
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where each ft g is the closure of ft e> 9 ft g is the boundary of ft g , ft^ is the finite ele- 
ment approximation of the closed bounded domain ft , and E is the total number of 
elements. 


Introducing the finite element discretization into equations (9) and (10) yields: 
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where 9ft e2 is the boundary of an element (ft g ) which lies on the part of the boundary 
(9ft2) on which the flux boundary condition is specified. 

The flux through inter-element boundaries should be continuous and the normal 
vectors at the interface of the two adjacent elements are in the opposite directions. 
Therefore all the inter-element fluxes in equation (12) cancel each other, and only 
the prescribed flux on the boundary 9ft 2 contributes to the final system of equations. 
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The consequence of eliminating the fluxes across inter-element boundaries in the finite 
element method is equivalent to enforcing the flux continuity condition across the 
inter-element boundaries . 

Inserting equations (3), (4), and (6) into equation (12) yields: 
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The global finite element system of equations are obtained by assembling the 
element system of equations. Therefore, the discrete finite element system of equa- 
tions are derived for an arbitrary element in the following discussions. 


Let <J> and <|> m be sets of basis polynomials to interpolate velocities and pres- 
sure respectively, i.e. , 


u. = u. <i> , sum on m, m = 1,M 

l im T m 


p = P n ^ n » sum on n, n = 1,N 


(15) 


where u. denotes the m-th nodal value of the velocity component u., p denotes the 

n-th nodal value of pressure, and M and N are the number of velocity nodes and the 
number of pressure nodes in an element, respectively. 

In the Galerkin finite element method, the test functions are selected from the 
same space of interpolating polynomials as the trial functions. Then and W can be 
expanded as : p 


W u “ * 

sum on k, k = 1,M 

W P = V* * 

sum on l, Z = 1,N 


where and b^ should vanish if the corresponding node lies on the part of the boundary 


on which the Dirichlet boundary condition is specified, respectively; otherwise, a^ and 
are arbitrary constants. 
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The finite element system of equations for an element is obtained by substituting 
equations (15) and (16) into equations (14) and (15), and is given as 
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where the subscript q ranges from 1 to M, and the arbitrariness of the coefficients a^ 
and b have been made use of in deriving equations (17) and (18). 


In matrix form, equations (17) and (18) are given as, for the three-dimensional 


case; 
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Uj is a column vector of nodal values of the velocity component u., p is a column 

vector of nodal pressure, <f> is a column vector of interpolating polynomials for 

velocity, \p is a column vector of interpolating polynomials for pressure, {b.c.} is a 

column vector contributed by the specified flux boundary condition, and the subscripts 
i and j range from one up to the number of spatial dimensions, respectively. 

For the two-dimensional case, the element system of equations can be obtained 
by deleting the appropriate third row and column sub-matrices from equations (19) 
and (20). 

The integrations in equations (21) through (25) were evaluated using the Gauss 
numerical quadrature method with three Gauss points for each coordinate direction. 

The assembled global system of equations was solved by a direct (Picard) iteration 
method using a frontal solver [1,12], and the solutions were updated using an under- 
relaxation method given as: 


Oj* = a a . n + (1 - a) a^ n ^ 


(26) 
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where a. represents any degree of freedom. 


a is the under- relaxation number, the 


superscripts n and n-1 denote iteration levels, 


and a.* is the updated solution. 


a 


(no under-relaxation) was used for all the flow variables for low Reynolds number 
flows; and a = 0.8 for velocities and a = 1 for pressure were used for high Reynolds 
number flows. 
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2-3. A New Mixed Interpolation for Velocity and Pressure 

The new flow elements used in the present finite element computation of Navier- 
Stokes equations are introduced in this section. For two-dimensional case, the 
velocities are interpolated using the biquadratic shape functions and the pressure is 
interpolated using the linear shape functions defined on a triangular element which is 
contained inside the quadratic element, as shown in Figure 1-a. The three pressure 
nodes are located at the three Gauss points of the three-point Gauss quadrature rule 
for quadrilateral elements [13], i.e., the same locations as those used in the Reduced 
Integration Penalty (RIP) method. The coordinates of the pressure nodes on the 
computational element are given as: 
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(-1//2 , 

-1//6 ) 

for n = 2 

(27) 


( 1//2 , 
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where £ = (£ ,n ) for 2-D case, and n denotes the pressure node numbers. The 

—n n n 

shape functions for each of the nodes are given as: 
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For the three-dimensional case, the velocities are interpolated by the tri- 
quadratic shape functions and the pressure is interpolated using the linear shape 
functions defined on a tetrahedral element which is contained inside the tri-quadratic 
brick element, as shown in Figure 1-b. The coordinates of the pressure nodes on 
the computational element are given as: 
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for n = 1 
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, 1//3 ) for n = 4 


( 29 ) 


where £ n = (S^t^.?^) for the 3-D case, and n denotes the node numbers of the 
pressure nodes. The shape functions for each of the pressure nodes are given as: 
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The shape functions given in equations (28) and (30) satisfy the relationship 

that : 




1 if £ = n 
0 if £ 4 n 


(31) 


where \p is the shape function for the £-th pressure node and C n denotes the coor- 
dinates of the n-th pressure node. Additional pressure interpolation polynomials tested 

T 

for two-dimensional flows were \p = (l,x,y). The latter interpolation scheme yielded 
the same computational results, up to several significant digits, as equation (28) for 
low Reynold number flows. As the Reynolds number is increased, i.e., Reynolds 
number of 10,000 case for the cavity flow, the pressure interpolation polynomials 
given in equation (28) were found to be the most stable numerically. 
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III. EXAMPLE PROBLEMS 


The finite element method described in the previous sections was tested and 
validated by solving a few example problems for which vast amounts of computational 
results and/or experimental data were available. These are: a lid driven cavity flow 

[5,14-18], a laminar backward- facing step flow [19-21], and a laminar flow in a square 
duct of strong curvature [22,23]. 

In the following discussions, solving the coupled momentum equation and the 
conservation of mass equation once is counted as an iteration. The convergence 
criterion used was 


a n+1 

|1 - "V" I < e , j = 1,T 
a. 

3 


(32) 


where a. denotes a nodal value of velocity component or pressure; T is the number of 

3 _4 

total degrees of freedom; and e is the convergence criterion, e = 1x10 was used for 

the two-dimensional flow cases; and e = 1x10 for the three-dimensional flow case. 

In elliptic flow computations, the velocity component may become vanishingly small 
near the stagnation regions, separation regions, and reattachment regions. Also 
pressure may become vanishingly small in some regions of the flow field. For these 

cases, i.e. , a^. n 0, it can be expected that equation (32) can hardly be satisfied. 

Therefore, equation (32) needs to be further supplemented by other constraint ( s) . 

For the present case, any velocity degree of freedom (a. n ) whose relative magnitude 

with respect to the maximum velocity in the flow domain is smaller than e^, i.e., 

|a. n /v I < e, and e, = 1x10 ^ , was screened out before equation (32) was applied, 
j riicix I l 

By the same way, any pressure degree of freedom (a. n ) whose relative magnitude with 

3 j| 

respect to the maximum pressure in the flow domain is less than i.e., aj /P max < 
e., was screened out from the convergence test. 


In the new mixed interpolation method introduced herein, the pressure is dis- 
continuous across element boundaries. Thus the nodal pressure at the velocity node 
was obtained by averaging all the pressure contributions made by the elements con- 
taining the node; and each of the contributions was computed using equation (28) for 
the two-dimensional case and equation (30) for the three-dimensional case. 


3-1. A Lid Driven Cavity Flow 

The cavity flow is described in Figure 2. No slip boundary condition, i.e., 
u = v = 0, was applied at all the boundaries except at y = 1 where u = 1 and v = 0. 
The Reynolds numbers considered were 400, 1000, 3200, 5000, 7500, and 10,000; where 
the Reynolds number is defined as R = pUL/y , U = 1 is the velocity of the lid, and 

L = 1 is the reference length. The Reynolds number was varied by adjusting the 
viscosity (y) of the fluid with the remainder of the parameters kept as constants. 

Two different grids, as shown in Figure 3, were used. 
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For the coarse grid case (Grid A), convergent solutions were obtained up to a 
Reynolds number of 3200. But the coarse grid could not resolve the smallest eddies 
(vortices 4 and 5 in Fig. 2), and thus the grid was not tested for Re > 3200. Velocity 
vectors obtained by using the coarse grid for Re < 1000 are shown in Figures 4-1 and 
4-b, respectively. The streamline contours for Re < 1000 obtained by using the coarse 
grid were qualitatively the same as those obtained by using the fine grid (Grid B), 
but the magnitude of the minimum stream function values at the center of the primary 
vortex were a few percent smaller than those obtained by using the fine grid. 

The velocity vectors (Re > 3200) and the streamline contours obtained by using 
the fine grid are shown in Figures 4 and 5, respectively. The streamline contour 
labels are given in Table 1. In Figure 5, it can be seen that the sizes of vortices 
4 and 5 of the present computational results are smaller than those of Schreiber and 
Keller [15] and Ghia et al. [16]. Details of the smallest eddies (vortices 4 and 5) 
need to be further studied using finer grid. 

The normalized pressure contours for all the Reynolds number cases obtained 
by using the fine grid are shown in Figure 6, and the pressure contour labels are 
given in Table 2. The normalized pressure (P) was obtained from the computed static 
pressure (p) using a relationship that P = pL ^/V ^/y » where L ref = 1 is the 

reference length, V » = 1 is the reference velocity, and y is the molecular viscosity 

of the fluid. re 

The horizontal velocity profiles at x = 0.5 are compared with various computa- 
tional results in Figure 7. For Re = 400, it can be seen that the minimum horizontal 
velocity due to Burggraf [ 14] is approximately 12 percent smaller than that of Ghia 
et al. [16], and that the same velocity due to Bercovier and Engelman [18] is approxi- 
mately 25 percent smaller than the same data. For Re = 1000, the horizontal velocity 
profile due to Bercovier and Engelman [16] further deviates from that of Ghia et al. 
or the present computational result. In general, the present computational results 
compare more favorably with those of Ghia et al. than those of Burggraf [14] and 
Bercovier and Engelman [16], as shown in Figure 7. 

For Re = 10,000, it can be seen that the horizontal velocity profile due to 

Schreiber and Keller [15] is slightly different from that of Ghia et al. [16], which 
shows that the global computational results can be slightly different depending on the 
order of difference approximation used for the boundary condition. The present 
computational results compare more favorably with those of Ghia et al. [16] than those 
of Schreiber and Keller [15]. But both of these finite difference solutions appeared 
only recently, and it is not quite clear which of these is better as yet. 

The lodal maximum or minimum values of the stream function at the center of 
vortices for the primary and the first three secondary vortices are compared with 
those of References 15, 16, and 17 in Table 3. It can be seen that the present 
computational results, obtained by using approximately one- fourth of the number of 
grid points used by others, compare favorably with the three data in general. 

Computation of the cavity flow was carried out in the order of increasing the 
Reynolds number. Uniform zero values for both velocity and pressure were used as 
an initial guess for Re = 10, the converged solution for Re = 10 was used as an 
initial guess for Re = 400, and so on. The required number of iterations were 5, 15, 
25, 42, 51, 60, and 194, for Reynolds numbers of 10, 400, 1000, 3200, 5000, 7500, 
and 10,000, respectively. The excessive number of iterations required for the Re = 
10,000 case, using the converged solution of Re = 7500 as initial guess, may be due 
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to the grid which is still not fine enough for the Reynolds number considered. The 
similar phenomenon was observed when the coarse grid (Grid A) was used to solve 
the Reynolds number of 3200 case. 

Study on the efficiency of the present numerical method has not been made as 
yet. But based on the required number of grid points and the required number of 
iterations to achieve comparable accuracy with those of the velocity -pres sure based 
finite difference flow solvers, the present method would not be unduly expensive 
compared with the finite difference methods. 


3-2. A Backward -Facing Step Flow 

A laminar backward- facing step flow is considered in the following. A descrip- 
tion of the problem is shown in Figure 8, and the experimental data for the flow can 
be found in Armaly et al. [19]. The Reynolds number, Re = pVD/y , is based on 
the hydraulic diameter (D = 0.0104 meters) and the bulk velocity (V = 0.6667 meters/ 
sec) at the inlet. The Reynolds numbers were varied by adjusting the molecular 
viscosity (y) of the fluid. For the Reynolds number less than approximately 450, 
there exists only one recirculation zone at the down- stream region of the backward- 
facing step. As the Reynolds number is increased beyond approximately 450, another 
recirculation zone appears at the top wall of the channel, the size of which grows 
further as the Reynolds number is increased. Experimental data show that a third 
recirculation zone appears at the bottom wall for the Reynolds number greater than 
approximately 1000. As the Reynolds number is increased beyond approximately 600, 
the three-dimensional effect becomes so strong that comparison between the two- 
dimensional computational result and the experimental data becomes less meaningful, 
which is discussed in detail in Armaly et al. [19]. Therefore, the computations were 
carried out from a Reynolds number of 100 up to a Reynolds number of 900 in the 
order of increasing the Reynolds number by 100 in the present study. For the coarse 
grid (Grid A) case, the uniform zero values for both velocity and pressure were used 
as an initial guess for the Re = 100 case, the converged solution of the Re = 100 case 
was used as an initial guess for the Re = 200 case, and so on. The required num- 
ber of iterations were 20, 40, 53, 89, 135, 161, 171, 172, and 183 for Reynolds 
numbers of 100 through 900, respectively. For the fine grid case (Grid B), the 
initial guess for all the Reynolds number cases were obtained by interpolating the 
coarse grid solutions using the multi- grid concept of Ghia et al. [16]. 

The computed velocity vectors and the streamline contours for Re = 400, 500, 
and 800 are shown in Figures 10 and 11. The streamline contour labels are given in 
Table 4. The onset of the second recirculation zone at the top wall can be seen in 
the streamline contour plot given in Figure 11-a, yet there exists no recirculation 
zone for Re = 400 as can be confirmed from the velocity vector plot in Figure 10-a. 

The computed static pressure was normalized using a relationship P = P L re f/ 

v re f/y> where L re f = 0.0049 meters is the step height and V re ^ is the bulk velocity 

at the inlet. The normalized pressure contours for Reynolds numbers of 400, 500, 
and 800 are shown in Figure 12. 

The size of recirculation zones versus Reynolds number is compared with 
experimental data in Figure 13. It can be seen that the computational results compare 
favorably with experimental data up to Reynolds number of approximately 600. 
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The top and bottom wall pressure for Reynolds numbers of 100 through 900 are 
shown in Figure 14. Neither the experimental data nor any other computational result 
for the wall pressure are available as yet, and no comparison could be made. 

Finite difference computation of the same backward-facing step flow can be 
found in Kim and Moin [20] and Kwak and Chang [21] among many others. The same 
level of agreement between the computational results and the experimental data as the 
present case can be found in Kim and Moin [20]. On the other hand, the streamline 
contour plot due to Kwak and Chang [21] did not show the top wall recirculation zone 
as clearly as that of the present case or that of Kim and Moin [20] . 


3-3. A Laminar Flow in a Square Duct of Strong Curvature 

A re-developing laminar flow in a square duct of strong curvature is considered 
below. The flow configuration is shown in Figure 15; and the experimental data for 
the flow is given in Humphery et al. [22], 

The Reynolds number based on the hydraulic diameter of 0.04 m and the bulk 
_ 2 

velocity of 1.98 x 10 m/sec was 790. Due to the symmetry of the flow domain, only 
one half of the duct was considered in the computation. The inlet plane was located 
2.8 hydraulic diameters upstream of the curved section, and the outlet plane was 8 
hydraulic diameters downstream. Discretization of the domain is shown in Figure 16. 
An analytical solution for the fully developed duct flow was prescribed at the inlet 
plane of the duct , and a vanishing normal stress boundary condition was prescribed 
at the exit. 

The computed velocity vectors at the three planes of the curved section are 
shown in Figure 17. It can be seen in Figure 17-c that there exists a weak recircula- 
tion zone extending up to 6 = 40° of the curved section. Computational results due 
to Humphery et al. [22] showed that the same recirculation zone extended beyond 
0 = 12°, and the computational results due to Rhie [23] showed that the recirculation 
zone extended beyond 6 = 30°. There does not exist fine grid computational results 
for the recirculation zone as yet. The secondary recirculation flows at three cross- 
sections are shown in Figure 18. The computational results showed that the grid used 
was not fine enough to resolve all the details of the flow field. It was also found that 
the frontal solver used in this study was not adequate to pursue further grid refine- 
ment. Fine grid computational results for the laminar flow in a square duct of strong 
curvature and other three-dimensional flow cases will be reported in the near future. 


IV. CONCLUSIONS AND DISCUSSION 


A new velocity -pressure integrated, mixed interpolation finite element method 
for the N avier- Stokes equations has been presented. The method has been tested 
for cavity flows with Reynolds number of 400 up to 10,000, a backward facing-step- 
flow of Armaly et al. [19], and a laminar flow in a square duct of strong curvature 
of Humpery et al. [22]. 

For the cavity flows, the present computational results compared favorably with 
those of Schreiber and Keller [15], Ghia et al. [16], and Gresho et al. [17], which 
were obtained by using approximately four times more grids than the present case. 
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For the backward- facing step flow, it was shown that the present computational 
method could capture the subtle pressure driven secondary recirculation zone at the 
top wall of the channel for Reynolds numbers greater than 500. The size of the 
recirculation zone also compared favorably with the experimental data. 

For the three-dimensional curved duct flow, grid refinement was not feasible 
due to the computer limitation and the frontal solver used in the present study. 
Detailed fine grid computation of the laminar flow in a square duct of strong curva- 
ture [22] will be reported in a near future together with other three-dimensional flow 
cases . 


It is usually known that the Bubnov- Galerkin method [24] (i.e., the test func- 
tions are selected from the same space as that of the trial functions, as in the present 
case) yields oscillatory solutions for high Reynolds number flows. An extensive dis- 
cussion on the finite element upwinding technique can be found in Brooks and Hughes 
[24] among many others. But the present computational results showed that no 
upwinding was necessary to obtain accurate solutions which are free of wiggles for 
the Navier- Stokes equations. 
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TABLE 1. STREAMLINE CONTOUR LABEL FOR CAVITY FLOW 


Label 

7 * 

Label 

i> 

Label 


A 

-0.11 

F 

-0.03 

K 

2.X10' 4 

B 

-0.10 

G 

-0.01 

L 

5.X10‘ 4 

C 

-0.09 

H 

-1.X10' 10 

M 

1.X10" 3 

D 

-0.07 

I 

1.X10' 6 

N 

2.X10' 3 

E 

-0.05 

J 

5.X10' 5 




V>: stream function 


TABLE 2. PRESSURE CONTOUR LABEL FOR CAVITY FLOW 


Label 



Reynolds 

Number (Re) 



400 

1000 

3200 

5000 

7500 

10000 

A 

-40. 

-100. 

-310. 

-470. 

-700. 

-900. 

B 

-30. 

-70. 

-220. 

-370. 

-500. 

-650. 

c 

-20. 

-50. 

-120. 

-270. 

-300. 

-400. 

D 

-10. 

-20. 

-70. 

-150. 

-150. 

-200. 

E 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

F 

10. 

30. 

100. 

120. 

300. 

400. 

G 

30. 

80. 

200. 

280. 

1000. 

1000. 

H 

50. 

— 

400. 

900. 


3000. 





TABLE 3. STREAM FUNCTION VALUES AT THE CENTER OF 
VORTICES FOR CAVITY FLOW 



Rey 

Schreiber* 
& Keller 

Ghia 0 
et. al. 

Gresho x 
et. al . 

present -1- 


400 

-0.1130 

-0.1139 

_ _ _ 

-0.1128 


1000 

-0.1160 

-0.1179 

-0.114 

-0.1169 


3200 

— 

-0.1204 

-0.118 

-0.1181 

rd <D 
E +-> 

5000 

— 

-0.1190 

-0.109 

-0.1173 

•1— S- 

s- o 

7500 

— 

-0.1200 

-0.108 

-0.1157 

Q_ > 

10000 

-0.1028 

-0.1197 

-0.101 

-0.1150 


400 

6.440xl0’ 4 

6 . 4235xl0' 4 



6 . 1810xl0’ 4 

r— 4 

1000 

1.700x10’ 3 

1 . 7510xl0’ 3 

1.760x10’ 3 

1 . 6594x10’ 3 

X 

3200 

— 

3 . 1396xl0’ 3 

3 . 290xl0’ 3 

2 . 6744x10’ 3 

Cl) 

+-> 

5000 

— 

3 . 0836xl0’ 3 

3 . 870x10’ 3 

2 . 7786x10’ 3 

S- 

o 

7500 

— 

3 . 2848xl0* 3 

4. 860xl0’ 3 

2 . 7396x10’ 3 

> 

10000 

2 . 960xl0’ 3 

3 . 4183x10" 3 

5. 540x10’ 3 

2 . 7528x10’ 3 


400 

1 . 450xl0’ 5 

1 . 4195xl0’ 5 

_ - - 

1 | 
1 . 3577x10’ 5 

CM 

1000 

2.170xl0’ 4 

2 . 3113xl0" 4 

2.0 xlO’ 4 

2 . 1951xl0’ 4 

X 

CL) 

3200 

1 

9 . 7820xl0’ 4 

1.20 xlO’ 3 

1 . 0465xl0’ 3 

+-> 

S- 

5000 

— 

1 . 3612xl0" 3 

1 . 490x10’ 3 

1 . 2675xl0’ 3 

o 

> 

7500 

— 

1 . 4671xl0" 3 

1 . 750x10’ 3 

1 . 3697x10’ 3 


10000 

— 

1 . 5183xl0" 3 

1 . 930x10’ 3 

1 . 4055x10’ 3 

CO 

3200 



7 . 2786xl0" 4 

5 . 860xl0’ 3 

6. 4440x10 ' 4 

i x 

<u 

5000 

— 

1.4564xl0’ 3 

1 . 230xl0’ 3 

1 . 3045x10’ 3 


7500 

— 

2 ,0462xl0" 3 

1 . 840xl0’ 3 

1 . 8426xl0’ 3 

i s 

10000 

— 

2.4210xl0" 3 

2 . 230x10’ 3 

2 . 1817xl0’ 3 


* 141x141 grids for Re - 400 and 1000, and 180x180 grids for 
Re=10 , 000 . 

o 257x257 grids for Re = 400, 5000, 7500, and 10,000, and 

129x129 grids for Re - 1000 and 3200. 

x 129x129 grids for Re<3200 and 257x257 grids for Re>5000. 

+ 65x65 grids (32x32 quadratic elements) for all the cases. 
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TABLE 4. STREAMLINE CONTOUR LABEL FOR 
BACKWARD-FACING STEP FLOW 


Label tp 

Label \p 

Label i> 

A -2.0X10' 4 

B -1 . 5X10' 4 

C -5 .0X10' 5 

D -1.0X1CT 5 

: E 0. 

F 1.0XKT 4 

G 5.0XI0' 4 

H 1.0X10' 3 | 

I 2. 0X10- 3 

J 3 . 0X10* 3 

K 3.467X10' 3 
L 3. 480X10- 3 
M 3. 50X10- 3 

1 


TABLE 5. PRESSURE CONTOUR LABEL FOR 
BACKWARD-FACING STEP FLOW 


Label 

Reynolds Number (Re) 

400. 

500. 

800. 

A 

-7.0 

-25.5 

-75.8 

B 

-6.0 

-23.0 

-75.0 

C 

-4.0 

-18.0 

- 73.5 

! D 

0.0 

-8.0 

-65.0 

E 

6.0 

0.0 

-45.0 

F 

12.0 

10.0 

-30.0 

; g 

20.0 

15.0 

-22.0 

; h 

28.0 

20.0 

-14.0 

! I 

32.0 

25.0 

-8.0 

i J 

33.5 

27.0 

0.0 

1 K ! 

34.5 

28.5 

3.0 


35.5 

29.5 

7.0 

i m 

36.0 

29.8 

10.0 

1 » 

36.5 

30.0 

11.5 

! 0 

36.8 

— 

12.0 

j 






Figure 2. Configuration, coordinates, and 
nomenclature of cavity flow. 









(b) Re =1,000, GRID A 



(c) Re = 3,200, GRIDB 


Figure 4. Velocity vectors for cavity flow. Re: Reynolds number 
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(d) Re = 5,000, GRID B 
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(e) Re = 7,500, GRID B 
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(f) Re = 10,000, GRID B 
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Figure 10. Velocity vectors for backward- facing step flow. 
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(a) Re = 400 x/h = 38 




(c) Re = 900 x/h = 38 


Figure 12. Pressure contours for the backward- facing step flow. 
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(b) BOTTOM WALL 


Figure 14. Wall pressure for backward- facing step flow 
fine grid (Grid B) solutions, Re = 100 - 900. 






Figure 16 . Discretization of the flow domain, (a) discretization of the 
cross-section, (b) grid in the flow direction. 
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Velocity vectors on the cur\ 
), (b) x /D 1/2 = 0.5, (c) x/ 
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Figure 18. Secondary recirculation flows, 
(a) 0 = 30°, (b) 9 = 60°, (c) 9 = 90°. 
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